## Script to run model and plots of rebel votes by votes and group
## Replicates Table 2 and Figure 1 in paper
## JBS last revised 3 August 2017

rm(list=ls())

library(plyr)
library(lme4)
library(stargazer)
library(arm)
library(mvtnorm)


##### SET YOUR PATH HERE ##############
##########################################

path <- "~/Dropbox (Personal)/Rebel Summaries/APSR_SKLLO_Repfiles/"

setwd(path)

load("CalcData/rebeldta9297.Rda")
load("CalcData/rebeldta0510.Rda")


mod1 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.9297[all.dta.9297 $Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod2 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.9297[all.dta.9297 $Party=="Con",],control=glmerControl(optimizer="bobyqa"))

mod3 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.0510[all.dta.0510 $Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod4 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.0510[all.dta.0510 $Party=="Con",],control=glmerControl(optimizer="bobyqa"))

summary(mod1)
summary(mod2)
summary(mod3)
summary(mod4)


groupsdta <- read.csv(paste(path,"RawData/groupmembership.csv",sep=""))

dta.w.groups9297 <- join(all.dta.9297,groupsdta,by="PubWhipID", type="left")
dta.w.groups9297$LCG[is.na(dta.w.groups9297$LCG)] <- 0
dta.w.groups9297$CORNER[is.na(dta.w.groups9297$CORNER)] <- 0
dta.w.groups9297$MONDAY[is.na(dta.w.groups9297$MONDAY)] <- 0
dta.w.groups9297$NTB[is.na(dta.w.groups9297$NTB)] <- 0
dta.w.groups9297$CONGROUPANY[is.na(dta.w.groups9297$CONGROUPANY)] <- 0
dta.w.groups9297$CONGROUPSUM[is.na(dta.w.groups9297$CONGROUPSUM)] <- 0

# Create one var for group membership to make tables work nicely with stargazer
dta.w.groups9297$group <- 0
dta.w.groups9297$group[dta.w.groups9297$CONGROUPANY == 1 | dta.w.groups9297$LCG == 1] <- 1
 

dta.w.groups0510 <- join(all.dta.0510,groupsdta,by="PubWhipID", type="left")
dta.w.groups0510$LCG[is.na(dta.w.groups0510$LCG)] <- 0
dta.w.groups0510$CORNER[is.na(dta.w.groups0510$CORNER)] <- 0
dta.w.groups0510$MONDAY[is.na(dta.w.groups0510$MONDAY)] <- 0
dta.w.groups0510$NTB[is.na(dta.w.groups0510$NTB)] <- 0
dta.w.groups0510$CONGROUPANY[is.na(dta.w.groups0510$CONGROUPANY)] <- 0
dta.w.groups0510$CONGROUPSUM[is.na(dta.w.groups0510$CONGROUPSUM)] <- 0

# Create one var for group membership to make tables work nicely with stargazer
dta.w.groups0510$group <- 0
dta.w.groups0510$group[dta.w.groups0510$CONGROUPANY == 1 | dta.w.groups0510$LCG == 1] <- 1


mod5 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups9297[dta.w.groups9297$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod6 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups9297[dta.w.groups9297$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

mod7 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups0510[dta.w.groups0510$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod8 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*group+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= dta.w.groups0510[dta.w.groups0510$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

summary(mod5)
summary(mod6)
summary(mod7)
summary(mod8)

#stargazer(mod1,mod2,mod3,mod4,type="latex",title="Effect of Government Status on Rebellion: Fractional Logit Models with Random Effects", align=T, column.labels=c("Labour","Conservative","Labour","Conservative"),dep.var.labels=c("1992--2001","2005--2010"),covariate.labels=c("In Government","25 Percent Most Rebellious","Majority","Leader","Gov*Rebellious"),style="ajps")

stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,type="latex",title="Effect of Government Status on Rebellion: Fractional Logit Models with Random Effects", align=T, column.labels=c("Labour","Conservative","Labour","Conservative","Labour","Conservative","Labour","Conservative"),style="ajps")


#### Create Plots of Simulations
#### Simulate backbenchers with mean victory margin (majority)

# Set number of simulations
sims<-1000

# Take random draws
dist<-rmvnorm(sims, fixef(mod5), as.matrix(vcov(mod5)))
dist2<-rmvnorm(sims, fixef(mod6), as.matrix(vcov(mod6)))
dist3<-rmvnorm(sims, fixef(mod7), as.matrix(vcov(mod7)))
dist4<-rmvnorm(sims, fixef(mod8), as.matrix(vcov(mod8)))


labdta9297 <- all.dta.9297[all.dta.9297$Party=="Lab",]
condta9297 <- all.dta.9297[all.dta.9297$Party=="Con",]
labdta0510 <- all.dta.0510[all.dta.0510$Party=="Lab",]
condta0510 <- all.dta.0510[all.dta.0510$Party=="Con",]

# Z is indicator of in gov or not
z<-c(0,1)

calc<-length(z)

lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
lo5<-numeric(calc)
hi5<-numeric(calc)
pe5<-numeric(calc)
lo6<-numeric(calc)
hi6<-numeric(calc)
pe6<-numeric(calc)
lo7<-numeric(calc)
hi7<-numeric(calc)
pe7<-numeric(calc)
lo8<-numeric(calc)
hi8<-numeric(calc)
pe8<-numeric(calc)

for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  pp5<-numeric(sims)
  pp6<-numeric(sims)
  pp7<-numeric(sims)
  pp8<-numeric(sims)
  
    for(j in 1:sims){
      #Labour (92-97); Moderate
      pp[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,7]*0*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T)+dist[j,5]*0+dist[j,6]*mean(labdta9297$start_year, na.rm=T)
      #Labour (92-97); extreme
      pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,7]*1*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T)+dist[j,5]*0+dist[j,6]*mean(labdta9297$start_year, na.rm=T)
      #Conservative (92-97); Moderate
      pp3[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*0+dist2[j,7]*0*z[i]+dist2[j,4]*mean(condta9297$Majority, na.rm=T)+dist2[j,5]*0+dist2[j,6]*mean(condta9297$start_year, na.rm=T)
      #Conservative (92-97); Extreme
      pp4[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*1+dist2[j,7]*1*z[i]+dist2[j,4]*mean(condta9297$Majority, na.rm=T)+dist2[j,5]*0+dist2[j,6]*mean(condta9297$start_year, na.rm=T)
      #Labour (05-10); Moderate
      pp5[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*0+dist3[j,7]*0*z[i]+dist3[j,4]*mean(labdta0510$Majority, na.rm=T)+dist3[j,5]*0+dist3[j,6]*mean(labdta0510$start_year, na.rm=T)
      #Labour (05-10); extreme
      pp6[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*1+dist3[j,7]*1*z[i]+dist3[j,4]*mean(labdta0510$Majority, na.rm=T)+dist3[j,5]*0+dist3[j,6]*mean(labdta0510$start_year, na.rm=T)
      #Conservative (05-10); Moderate
      pp7[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*0+dist4[j,7]*0*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T)+dist4[j,5]*0+dist4[j,6]*mean(condta0510$start_year, na.rm=T)
      #Conservative (05-10); Extreme
      pp8[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*1+dist4[j,7]*1*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T)+dist4[j,5]*0+dist4[j,6]*mean(condta0510$start_year, na.rm=T)
    }
    
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  pe5[i]<-invlogit(mean(pp5))
  lo5[i]<-invlogit(quantile(pp5, 0.025))
  hi5[i]<-invlogit(quantile(pp5, 0.975))
  
  pe6[i]<-invlogit(mean(pp6))
  lo6[i]<-invlogit(quantile(pp6, 0.025))
  hi6[i]<-invlogit(quantile(pp6, 0.975))
  
  pe7[i]<-invlogit(mean(pp7))
  lo7[i]<-invlogit(quantile(pp7, 0.025))
  hi7[i]<-invlogit(quantile(pp7, 0.975))
  
  pe8[i]<-invlogit(mean(pp8))
  lo8[i]<-invlogit(quantile(pp8, 0.025))
  hi8[i]<-invlogit(quantile(pp8, 0.975))
  
  print(i)
}



#92-2001
pdf("FinalPlots/Reb_by_Reb_92_2001.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe1, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.065), ylab="Probability of Rebellion", xaxt="n", col="grey80", xlab="")
segments(0.15, lo1[1], 0.15, hi1[1], lwd=2, col="grey80")
segments(0.65, lo1[2], 0.65, hi1[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe2, pch=16, cex=1.5, col="grey80")
segments(0.20, lo2[1], 0.20, hi2[1], lwd=2, col="grey80")
segments(0.70, lo2[2], 0.70, hi2[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe3, pch=15, cex=1.5, col="grey20")
segments(0.25, lo3[1], 0.25, hi3[1], lwd=2, col="grey20")
segments(0.75, lo3[2], 0.75, hi3[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe4, pch=16, cex=1.75, col="grey20")
segments(0.30, lo4[1], 0.30, hi4[1], lwd=2, col="grey20")
segments(0.80, lo4[2], 0.80, hi4[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Non-Group MPs", "Ideological Group", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()

#2005-2010
pdf("FinalPlots/Reb_by_Reb_01_2005.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe5, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.065), ylab="Probability of Rebellion", xaxt="n", col="grey80", xlab="")
segments(0.15, lo5[1], 0.15, hi5[1], lwd=2, col="grey80")
segments(0.65, lo5[2], 0.65, hi5[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe6, pch=16, cex=1.5, col="grey80")
segments(0.20, lo6[1], 0.20, hi6[1], lwd=2, col="grey80")
segments(0.70, lo6[2], 0.70, hi6[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe7, pch=15, cex=1.5, col="grey20")
segments(0.25, lo7[1], 0.25, hi7[1], lwd=2, col="grey20")
segments(0.75, lo7[2], 0.75, hi7[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe8, pch=16, cex=1.75, col="grey20")
segments(0.30, lo8[1], 0.30, hi8[1], lwd=2, col="grey20")
segments(0.80, lo8[2], 0.80, hi8[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Non-Group MPs", "Ideological Group", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()

